Stability of Metal Nanowires at Ultrahigh Current Densities 
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We deveiop a generaiized grand canonicai potentiai for the baiiistic nonequiiibrium eiectron dis- 
tribution in a metai nanowire with a finite appiied bias voitage. Couiomb interactions are treated in 
the seif-consistent Hartree approximation, in order to ensure gauge invariance. Using this formalism, 
we investigate the stability and cohesive properties of metallic nanocylinders at ultrahigh current 
densities. A linear stability analysis shows that metal nanowires with certain magic conductance 
values can support current densities up to 10 11 A/cm 2 , which would vaporize a macroscopic piece 
of metal. This finding is consistent with experimental studies of gold nanowires. Interestingly, our 
analysis also reveals the existence of reentrant stability zones — geometries that are stable only under 
an applied bias. 



PACS numbers: 61.46.+W, 68. 65. La 47.20.Dr, 66.30.Qa 



I. INTRODUCTION 

Metal nanowires have been the subject of many ex- 
perimental and theoretical studies, both for their unique 
properties and potential applications (see Ref. [lj for a 
review of the field). One of the most remarkable prop- 
erties of metal nanowires is their ability to support ex- 
tremely high current densities without breaking apart or 
vaporizingi2»MiMiL2i For noble metals, experiments can 
be carried out in air, and the first few peaks in conduc- 
tance histograms can withstand applied voltages as high 
as one volt, and even two volts for the first peak, cor- 
responding to one conductance quantum Go = 2e 2 /h. 
Let us estimate the corresponding current density: For 
a ballistic metallic conductor in the form of a cylinder 
of radius R, the electrical conductance G is given ap- 
proximately by the Sharvin formula G ~ Go(kpR/2) 2 , 
where kp is the Fermi wavevector. Therefore the current 
density at applied voltage V is 
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where n is the number density of conduction electrons, 
vf is the Fermi velocity, and sf is the Fermi energy. For 
an applied bias of a few Volts, the factor eV7 ef is of or- 
der unity, and the current density is of order 10 n A/cm 2 . 
Such high current densities would vaporize a macroscopic 
wire, thus prompting questions on the reason for the re- 
markable stability of metal nanowires. 

The first part of the answer to this question is that 
metal nanowires are typically shorter than the mean-free 
path Lj n for inelastic scattering, so that the conduction 
electrons can propagate through the wire without gen- 
erating excitations such as phonons^ that heat the wire. 
Instead, most of the dissipation takes place in the macro- 
scopic contacts for the outgoing electrons. However, 
the absence of equilibration of the electron distribution 
within the nanowire raises another, more fundamental, 
question: What is the effect of a highly nonequiiibrium 
electron distribution on the stability of a metal nanowire, 
given that the conduction electrons play a dominant role 



in the cohesion of metals? That is the question to which 
the present article is devoted. 

Under a finite bias, the scattering states of right- and 
left-moving electrons in a nanowire are populated differ- 
ently, even if there is no inelastic scattering within the 
wire. An adequate treatment of the electron-electron in- 
teractions is crucial to correctly describe this nonequiiib- 
rium electron distribution. Some studies of transporwas* 
and cohesionii in metal nanowires at finite bias did not 
include electron-electron interactions, so that the calcu- 
lated transport and energetics depended separately on 
both the left and right chemical potentials and 
thus violating the gauge invariance condition: The cal- 
culated physical quantities should depend only on the 
voltage difference eV — \i+ ~ and should be invari- 
ant under a global shift of the electrochemical potential, 
since the total charge is conserved^ A self-consistent 
formulation of transport and cohesion at finite bias has 
recently been developed based on ah initio and tight- 
binding methods These computational tech- 
niques are particularly well-suited to the study of atomic 
chains, but can become intractable for larger nanostruc- 
tures. An analytical approach to this problem is needed 
to study the interesting mesoscopic effect a 17 . 18 which oc- 
cur in systems intermediate in size between the macro- 
scopic and the atomic scale. 

In this paper, we extend our continuum 
model 18 i 19 . 20 i 21 . 22 i 23 . 24 i 25 i 26 of metal nanowires to 
treat the ballistic nonequiiibrium electron distribution 
at finite bias. Our model provides a generic description 
of nanostructures formed of simple, monovalent metals. 
It is especially suitable for alkali metals, but is also 
appropriate to describe quantum shell effects due to 
the conduction-band s electrons in noble metals. For 
a fuller discussion of the domain of applicability of 
our continuum approach, see Ref. l22t In the present 
work, Coulomb interactions are included in the self- 
consistent Hartree approximation, in order to ensure 
gauge invariance. 

For a system out of equilibrium, there is no general way 
to define a thermodynamic free energy. By assuming that 
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the electron motion is ballistic, however, the energetics of 
the biased system can still be described by a nonequilib- 
rium free energy^ which can be used to study the stabil- 
ity and cohesion of nanowircs at finite bias. We find that 
metal nanocylinders with certain magic conductance val- 
ues, G/Gq = 1, 3, 6, 12, 17, 23, 34, 42, . . ., can support cur- 
rent densities up to 10 11 A/cm 2 . Our finding is consistent 
with experimental results for gold nanocontacts^iiSiSil 
(G < 5Go) and atomic chains^ (G ~ Go), but implies 
that the magic wires with G > 5Gq are also extremely 
robust. Furthermore, we predict a number of nanowire 
geometries that are stable only under an applied bias. 

This paper is organized as follows: In Sec. |Hj we de- 
velop a formalism to describe the nonequilibrium thermo- 
dynamics of a mesoscopic conductor at finite bias. In Sec. 
IIIII we apply this formalism to quasi-one-dimensional 
conductors, and obtain gauge-invariant results for the 
Hartree potential, grand canonical potential, and cohe- 
sive force of metal nanocylinders at finite bias. In Sec. II VI 
we perform a linear stability analysis of metal nanocylin- 
ders at finite bias, the principal result of the paper. Sec- 
tion presents some discussion and conclusions. Details 
of the stability calculation are presented in Appendix 1X1 



II. SCATTERING APPROACH TO 
NONEQUILIBRIUM THERMODYNAMICS 

We consider a metallic mesoscopic conductor con- 
nected to two reservoirs at common temperature T = 
with respective electrochemical potentials 
fj,± = ef + eV±, where ep is the chemical potential for 
electrons in the reservoirs at equilibrium, e is the electron 
charge, and V+(_j is the voltage at the left (right) reser- 
voir. Because the screening of electric fields in metal 
nanowires with G > Gq is quite good, the presence of 
additional nearby conductors (such as a ground plane) 
has a negligible effect on the transport and energetics of 
the system, and is therefore not considered. 

While there is no general prescription for constructing 
a free energy for such a system out of equilibrium, it is 
possible to do so based on scattering theory 2 ^ if inelastic 
scattering can be neglected, i.e., if the length L of the con- 
ductor satisfies L <C L in . In that case, scattering states 
within the conductor populated by the left (right) reser- 
voir form a subsystem in equilibrium with that reservoir. 
Dissipation only takes place for the outgoing electrons 
within the reservoir where they are absorbed. Treating 
electron-electron interactions in mean-field theory, it is 
then possible to define a nonequilibrium grand canonical 
potential O of the system, 

n[n + , f i-,u(f)} = n [ t i + ,fx-,u(n} 
i 



d 6 r[p-(n+P+(r)]U{f), (2) 



where £Iq describes independent electrons moving in the 
mean field U, p± are the number densities of electrons 



(— ) and of ionic background charges (+), and the second 
term on the r.h.s. of Eq. (J2J corrects for double-counting 
of interactions in Slo- Since the electrons injected from 
the left and right reservoirs are independent, aside from 
their interaction with the mean field, Qo is given by the 
sum 



Slo\p+,(X-,U(f)] 



a=± 



^qK, U(f)], 



where 



n a \Ma,U] = ~k B T dE g a {E)\n 1 + e 



-p{E-fi a ) 



(3) 



(4) 



is the grand canonical potential of independent electrons 
moving in the potential U, in equilibrium with reservoir 
a, and the injectivity&ZS^ 



-w-3sl>{^-^.} 



(5) 



is the partial density of states of electrons injected by 
reservoir a. Here S ia = S 7a [E,U(f)} is the submatrix 
of the electronic scattering matrix describing electrons 
injected from reservoir a and absorbed by reservoir 7, 
and is a functional of the mean-field potential. 

The number density p~{r) of the conduction electrons 

is 



P-(r) = 



5Q 
SU(f) 



= fdEg a (r,E)f(E,n a ), (6) 

a=± J 



where f(E,p) = (1 + exp[f3(E — p,)]) is the Fermi- 
Dirac distribution function, and 



la (f ' E) - ~ hi ^ Tr { sl < a W(F) - sufj ^ } 



(7) 



is the local partial density of state^^^ for electrons 
injected from reservoir a. In Eqs. © and 0, S/6U(r) 
denotes the functional derivative. 

The mean-field potential U is determined in the 
Hartree approximation by 



U(f) = d 3 r'V{r-f")5p(r'), 



(8) 



where Sp(f) = p_ (r) — p+{r) is the local charge imbal- 
ance in the conductor and V(r ) = e 2 /\f\ is the Coulomb 
potential. The Hartree potential depends on the electro- 
chemical potentials of the left and right reservoirs. 

The whole formalism |PJ Q is very similar for any 
mean-field potential that is a local functional of the elec- 
tron density^ but we choose to work with the Hartree 
potential for simplicity. The exchange and correlation 
contributions to the mean field are taken into account 
in the present analysis only macroscopically, 22 by fixing 
the background density p + to its bulk value. Throughout 
this paper, we assume p + = k F /3n 2 — const, within the 
conductor (jellium model). 
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III. QUASI-ONE-DIMENSIONAL LIMIT 

Equations J3J JHJ provide a set of equations at finite 
bias that have to be solved self-consistently. For a con- 
ductor of arbitrary shape, these equations may be quite 
difficult to solve. We therefore restrict our consideration 
in the following to quasi-one-dimensional nanoconduc- 
tors, with axial symmetry about the z-axis. The shape of 
the conductor is specified by its radius R(z) as a function 
of z, and we assume R(z) <C L. For such a quasi-one- 
dimensional geometry we can approximately integrate 
out the transverse coordinates, replacing the Coulomb 
potential by an effective one-dimensional potential 



V(z,z') = 



\(z-z') 2 + \ (R 2 (z) + R^iz'))} 1 / 2 ' 



(9) 



where r\ is a parameter of order unity. The longitudinal 
potential V(z, z 1 ) has to be supplemented with a trans- 
verse confinement potential, which we take as a hard wall 
at the surface of the wire^i This boundary condition ne- 
cessitates a careful treatment of surface charges, as dis- 
cussed below. As a consistency check, our final results for 
the stability and cohesion are independent of the value 
of r\ chosen in the effective Coulomb potential. 

With this form of the Coulomb potential, the mean 
field U(z) becomes a function of the longitudinal coordi- 
nate only, and Eqs. ©~© reduce to a series of one- 
dimensional integral equations, which are much more 
tractable. The grand canonical potential of a quasi- 
cylindrical wire of length L is 

fi[{/x±}, R(z), U{z)] = fl [{p ± }, U(z)} 



dz[p.(z)+p + (z)]U(z), (10) 



2 Jo 



where fio is still given by Eqs. (JSJ-JSJ, with Q a 
Q a [p Ja ,R(z),U(z)]. Here 



sn 



JdEg a (z,E)f(E,p a ) (11) 



is the linear density of conduction electrons, where 
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ya. 



SU(z) SU(z) 



is the injectivity of a circular slice of the conductor at 
z. The scattering matrix S — S[E, R(z), U(z)] is now 
a functional of R(z) and U(z). In order to compensate 
for the depletion of surface electrons due to the hard- 
wall boundary condition, the linear density of positive 
background charges is taken to be 



P+( z ) 



k 3 F R{z) 2 k 2 F R{z)^/l + R 2 (z) 



3tt 



3tt 



R{z)R zz {z) 



(13) 



where R z {z) = dR{z)/dz and R zz {z) = d 2 R(z)/dz 2 . 
The second term on the r.h.s. of Eq. (|13|) corre- 
sponds to the well-known surface correction in the free- 
electron models The last term represents an integrated- 
curvature contribution, which is found to be a small cor- 
rection. The prescription given in Eq. 1)13(1 is essentially 
equivalent to the widely employed practice of placing the 
hard-wall boundary at a distance d = 37r/8fcp outside 
the surface of the metali^i 
The Hartree potential is 



U(z) 



dz V(z, z')Sp(z'), 



(14) 



where Sp(z) — p~(z) — p+(z). 

Equations l|§|l- 114(l provide a natural, gauge- 
invariant, generalization of the nanoscale free-electron 
model, 1 ® which has been successful in describing many 
cquilibriumiii^ and linear-responsei&i&S! properties of 
simple metal nanowires, to the case of nanowires at finite 
bias. This formalism represents a considerable simpli- 
fication compared to ab initio approachesii^i^i^i^ or 
even traditional jellium calculationSf^S and permits 
analytical results for the cohesion and stability of metal 
nanocylinders at finite bias. 



Quasicylindrical nanowire without backscattering 



Consider a nearly cylindrical nanowire, with radius 



R(z) = R + XSR(z), 



(15) 



where A is a small parameter and 5R(z) is a slowly- 
varying function. The couplings of the nanowire to 
the reservoirs are assumed ideal, so that electrons en- 
ter or exit the conductor without backscattering. For 
sufficiently small A, electron waves partially reflected or 
transmitted by the small surface modulation can be ne- 
glected, because they give a negligible contribution to 
the density of states. The right (left)-moving electrons in 
the bulk of the wire are thus in equilibrium with the left 
(right) reservoir. Moreover, the Hartree potential U{z) 
varies slowly with z, and can be taken as a shift of the 
conduction-band bottom in the adiabatic approximation. 
The injectivity therefore simplifies to 

g a (z,E) = ^g(z,E-U(z)), 

where g(z, E) is the local density of states for free elec- 
trons in a circular slice of radius R(z). Eq. ||2J can be 
rewritten as 



n a [p,R(z),U(z)} = 1 I dz 



H-U{z) 



dE 



x[E- p + U(z)]g T (z,E), (16) 

where the integration variable E is no longer the total 
energy of an electron, but rather its kinetic energy. Here 
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g T (z, E) = -J dE' g(z, E')df(E - E', 0)/dE is a convo- 
luted local density of states in a slice of the wire, and 
can be used to obtain finite-temperature thermodynamic 
quantities from their zero-temperature expressions^ so 
that Eq. (|16() is equivalent to the usual definition of the 
grand canonical potential Q). Similarly, the linear den- 
sity of electrons can be written in terms of the convoluted 
density of states as 

M*) = ^ = ~£/ a U(Z) dE gT (z,E). (17) 
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The convoluted density of states of a circular slice 
of the nanowire can be expressed semiclassically as^ 
gr(z, E) = §t{z, E) + 5gp{z 1 E), where gp is a smoothly- 
varying function of the geometry, known as the Weyl 
term, and 5gp{z,E) is an oscillatory quantum correc- 
tion. The temperature dependence of the Weyl term is 
negligible, g T = g x (l + 0(T/T F ) 2 ) , where T F is the 
Fermi-temperature. The zero-temperature value is 

_, kldV(z) [E 

9 ^ E) = ^p-\j- F _ 

fcjggCg) , k F 8C(z) fef 
87T£ F + 6tt 2 £ F V E ' 1 j 

where 9V(z), dS(z) and dC(z) are, respectively, the vol- 
ume, external surface-area, and external mean-curvature 
of a slice of the wire. The fluctuating part Sgp can be 
obtained through the trace formulaSi 

k 2 x - 

8g T {z, E) = 22 fv W a vw {T) 

WV 

x^4^co S vw (z,E), (19) 

where the sum is over all classical periodic orbits (v,w) 
in a disk billiardSSiiSi of radius R{z). Here the factor 
= 1 for v = 2w, and 2 otherwise, accounts for 
the invariance under time-reversal symmetry of some or- 
bits, L vw {z) = 2vR(z)sm(irw/v) is the length of peri- 
odic orbit (v,w), 9 vw (z,E) — k F L vw (z)^J E j e F — 3i>7r/2, 
and a vw (T) = t vw / smlir^ is a temperature-dependent 
damping factor, with t vw = nk F L vw T/2T F . 

The semiclassical approximation, Eqs. (|TH|) and (|TU|l . 
allows for an analytical solution for the ballistic nonequi- 
librium electron distribution in a metal nanowire at fi- 
nite bias. It also enables us to carry out a linear stability 
analysis of metal nanowires at finite bias, with analytical 
results for the stability coefficients. Although these cal- 
culations could in principle be carried out using a fully 
quantum mechanical solution of the electronic scatter- 
ing problem, the semiclassical approximation has been 
shownSiiffl to accurately describe the long-wavelength 
surface perturbations that are the limiting factor— in 
the stability of long nanowires. 



FIG. 1: (Color online) The Hartree potential Uq for electrons 
in a cylindrical nanowire at finite temperature T = 0.008 Tf, 
with a symmetric potential drop, versus the radius Ro of the 
wire. 

Eqs. ifTil) and (|T7J) provide a set of self-consistent equa- 
tions to solve for the ballistic nonequilibrium electron dis- 
tribution in a quasi-one-dimensional nanoconductor at fi- 
nite bias. Once the distribution p_ is obtained, the grand 
canonical potential fl of the electron gas may be calcu- 
lated from Eqs. (|10f> and (|16fl . The functional depen- 
dence of il[R(z)] yields information on the cohesioniSiSl 
and stabilitj*2ii22i24i£ of a metal nanowire, as in the equi- 
librium case. 

Solution for a cylindrical nanowire; Hartree 
potential and tensile force 

For an unperturbed cylinder, the mesoscopic Hartree 
potential Uq that simultaneously solves Eqs. 114|) and 
(|T7jl is only a function of the radius Rq, voltage eV = 
/(+ — fi-, and temperature T, and is constant along the 
wire, neglecting boundary effects, which are important 
only within a screening length ( ~ fc^ 1 ) of each con- 
tact. (This description is valid for wires with 1 <C kph < 
k F L- m .) Uq is independent of the choice of r\ in the 
Coulomb interaction, Eq. @, and can be determined by 
the charge neutrality condition 

Q = ~e[JV_0i+ - Uq) + N_(ii_ - U Q )] - eN+ = 0, (20) 

where \ N- (/i = fi± — U ) = \ J Q L dz J^dE gr(z, E) is the 
number of right(left)-moving electrons in the cylindrical 
wire, and JV+ is the total number of background positive 
charges. Equation (12011 gives a relation^ 

U Q = U^ s \Rq,V,T) + i(/i+ + m-) - e F , (21) 

(s) 

where Uq is calculated with a symmetric voltage drop 
V+ = —V- = \V . Equation (|2TJl guarantees that all 
physical properties of the system calculated in the fol- 
lowing are just functions of the voltage V, and not of /x+ 
and /!_ separately. 
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Using these expressions, one can solve Eq. (|2T))) for Uq 
for a symmetric potential drop, /i± = Ep ± \eV ■ The 
solution is shown in Fig.^ as a function of radius Ro at 
two different voltages V. In the equilibrium case (V = 0), 
/To oscillates about zero, exhibiting cusps at the subband 
thresholds, and increasing in amplitude as Ro decreases, 
due to the quantum confinement. Note that in equilib- 
rium, Uq — > as Rq — > oo, consistent with the well-known 
behavior of bulk jellium. At finite bias, each cusp in Uq 
splits in two, corresponding to the subband thresholds for 
left- and right-moving electrons: 



(i± = e„(R ) + U (R ,V,T), 



(22) 



where £„(i?o) are the eigenenergies of a disk billiard of 
radius Rq. This is illustrated in Fig. ^ for eV = 0.5ef- 
Note that in addition to the splitting, there is a substan- 
tial shift of the peak structure at finite bias. 

Using Eq. i|10|) . the grand canonical potential of a 
cylinder is now found to be 



n(R , V, T) = n [{n±}, R , U ] - N+U Q . 



(23) 



Note that is invariant under a global shift of the poten- 
tial Uq, due to an exact cancellation in the two terms on 
the r.h.s. of Eq. i|23|) . The tensile force in the nanowire 
provides direct information about cohesion, and is given 
by 



F(Ro,V,T) 



on 



RlL,V.T 



Figure [21 shows the tensile force of a metal nanocylindcr 
as a function of its cross section for two different bias 
voltages. To facilitate comparison with the stability di- 
agrams in Sec. IIVI below, the cross section is plotted in 
terms of the corrected Sharvin conductance^ 



G s = G Q [(k F R Q /2) 2 - k F R /2], 



(24) 



which gives a semiclassical approximation to the elec- 
trical conductance. In Fig. F < corresponds to 
tension, while F > corresponds to compression. As 
shown below, the cusps in the cohesive force at the sub- 
band thresholds correspond to structural instabilities of 
the system. 

In Fig. the force calculated at zero bias is very sim- 
ilar to previous results^SdL^ based on the free-electron 
model, even though those calculations did not respect the 
charge neutrality (|2C)|> enforced by Coulomb interactions. 
The reason for the good agreement is that the contribu- 
tion of the Hartree potential to the energy of the sys- 
tem is a second-order mesoscopic effect at zero biasj22i2£ 
which is essentially negligible for G > 3Gq. An earlier 
calculation^ that did invoke charge neutrality obtained 
a very different — and incorrect — result, because the sec- 
ond term on the r.h.s. of Eq. (|23l) was omitted, resulting 
in a double-counting of Coulomb interactions. 

Figure [21 shows that the cohesive force of a metal 
nanowire can be modulated by several nano-Newtons for 
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eV=0.5s w 




4 5 6 
G s [2e*/h] 

FIG. 2: (Color online) Tensile force in a metal nanocylinder 
versus cross-sectional area, at two different bias voltages at 
temperature T = 0.008 Tf- The cross section is plotted in 
terms of the Sharvin conductance Gs, Eq. JSJ, and the force 
is given in units of ef/^f (Af being the Fermi wavelength), 
which is 1.7 nN for Au. 



a bias of a few Volts. Such a large effect should be observ- 
able experimentally using appropriate cantilevers^^^ 
although the intrinsic behavior might be masked by elec- 
trostatic forces in the external circuit. In contrast to 
the case at V = 0, the Coulomb interactions play an 
essential role in determining the cohesive force at finite 
bias, since the positions of the peaks depend sensitively 
on the Hartree potential Uo of the ballistic nonequilib- 
rium electron distribution, shown in Fig. ^ The gauge- 
invariant result shown in Fig.[2thus differs substantially 
from previous results jii where screening was not treated 
self-consistently. Gauge-invariant results for the nonlin- 
ear transport through metal nanocylinders will be pre- 
sented elsewhere 4£ 



IV. LINEAR STABILITY OF A CYLINDER AT 
FINITE VOLTAGE 

In this section, we perform a linear stability 
analysis2ii22 for cylindrical nanowires under a finite bias 
V. Coulomb interactions are included self-consistently 
using the formalism of Sec. II I II Although the details 
of the calculation are rather complicated, the method 
is conceptually straightforward: Having found the self- 
consistent solution (|2*3^l for a cylinder, we perturb the 
cylinder as in Eq. (|15|) . and expand the free energy up 
to second order in the small parameter A. First, the 
self-consistent integral equations l|14|) and i|17|) for the 
ballistic noncquilibrium electron distribution are solved 
using first-order perturbation theory in A. Then the free 
energy is calculated using Eqs. (|10|l and (|16|) . 

The radius of the wire is given by Eq. (|15fl . with a 
perturbation function 



5R(z) = J2b(Q> iqz 



b(q)* = b(- q ). 



The surface perturbation XSR(z) is subject to a 
constraint fixing the total number of atoms in the 
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wire. In previous works, we have considered various 
constraints )22i2i which allow one to adjust the surface 
properties of the wire to model various materials. The 
simplest such constraint is volume conservation, under 
which the coefficient 6(0) is fixed by 

A 



Other reasonable constraints do not lead to qualitatively 
different conclusions. 

Within linear response theory, we can expand Sp(z) 
around Uo to linear order in U(z) — Uq, where Uq is the 
mesoscopic Hartree potential for the corresponding un- 
perturbed cylindrical wire. One gets 



5p(z) ~ 5po(z) ~ J dz 
dz\dz 2 



, Sp-(z) 



SU(z') 
Sp-{z 



SU{ Zl ) 



U 



V(z 1 ,z 2 )Sp(z 2 ), (25) 



where Eq. I|14|l has been used and Spo(z), called the bare 
charge imbalance, is defined as 

S Pa(z)=~J2 / dEg T (z,E)-p + (z). (26) 

Now defining the dielectric function 

5p-{zi) 



e(zi,z 2 ) = 5(z 1 -z 2 )- 



dz 3 



Uo 



SU(z 3 ) 

we can rewrite Eq. H25J) as 

Sp(z) = dz'e~ 1 (z, z')Spo(z'), 



V(z 3 ,z 2 ), (27) 



(28) 



where e 1 is the inverse dielectric function which satisfies 
/ dz 3 e~ 1 (zi,z 3 )e(z3, z 2 ) = S(zi - z 2 ), and 



Spo(z) = 6p {z) - / dzi 



Sp-(z) 



5U(zi) 



Uo 



(29) 



U 



The functional derivative 5p-(z)/SU(z') can be calcu- 
lated using Eq. (|I7|I . and is found to be 



dU(z') 



S(z - z') 



U 



^2 gr(z,p a - U ). 



Now we can expand the nonequilibrium grand canoni- 
cal potential (|I0|l as a series in A. In order to do so, we 
first expand Eq. (^DJ) around Uq. Using Eqs. ifl^jl. (fTTjl . 
d2Hl> and Jini), one gets 



n = n [{p±},R(z),U ] -Uo / dz5p (z) 

Jo 

1 f L _ 

+ - dz dz' Spo(z)V(z, z')5po(z') 

2 Jo 

+ ^T~J2 f L dzg T (z,e a )-N + U , (30) 

a = ± J° 



where e± = p± — Uq, and the screened potential V(z, z') 
is defined as 



V(z,z')= / dziV{z,Zi)e- x {z u z'). 



(31) 



At this point, all quantities have been expressed in 
terms of the local density of states gr{z,E) and the 
Coulomb interaction V(z, z'), whose expansions in series 
of A are presented in Appendix [S] In the end, the ex- 
pansion of the grand canonical potential as a series in A 
is found to be 

fi = fl(i? , V, T) + \ 2 L S(g; Ro, V, T)\b(q)\ 2 (32) 

plus terms C(A 3 , L°), where f2(i?o, V, T) is given by Eq. 
(|23|l . and the mode stiffness 



E(q;Ro,V,T) = a(q;Ro,V,T) 
+ Rc 



^V>-^V) 2 . 03) 



Here V(g) and e(q) are, respectively, the Fourier trans- 
forms of V^(z) and e^(z), the Coulomb potential and 
dielectric function of an unperturbed cylinder. The fac- 
tors g± and g^ are given in Eqs. i|A6(l and l|A7() . 
The factor a(q; Ro, V, T) comes from the expansion of 
^o[{m±}, R{z), Uo] and is found to be 



a{q;Ro,V,T) 



ir<7 s e\ + e 2 _ 
Ro e 2 F 



-KGsRo- 



+ 



( d 2 



1 d 



\9Rq Ro ORq 



3/2 . 3/2' 

r + + e_ 
372 



V sheU (R ,V,T), (34) 
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FIG. 3: (Color online) Total mode stiffness E(q = 0;R , V, T), 
and leading-order contribution a(q = 0; Ro, V, T), at zero and 
finite bias, for T = 0.008 Tf. 
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FIG. 4: (Color online) The factor D(Rq, V), coefficient of the 
q 2 contribution to H(g; Ro, V), at temperature T = 0.008 Tp. 
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where cr s is the surface tension, 7 S is the curvature en- 
ergy, and Kheii(-Ro) V, T) is the mesoscopic electron- shell 
potential, 23 given self-consistently at finite bias by 



V shen (R ,V,T) = — 

a — ± wv 



x 4 C0S (k%L vw - 3vtt/2) , (35) 



where kp = kpy/ e±/eF- In the present article, the 
values^ a s = ephp/Wir and 7 S = Ief^f /9tt 2 , appropri- 
ate for a constant- volume constraint, are used through- 
out. Inserting material-specific values^ does not lead to 
a significant change in the stability diagram. Equations 
(1321 l|35|) represent the central result of this paper. 

Figure |2| shows the long-wavelength mode stiffness 
S(g = 0) at zero and finite bias. For comparison, the 
leading-order contribution a(q = 0) is plotted as a dashed 
curve. The second term on the r.h.s. of Eq. (|33[) . which is 
second-order in the induced charge imbalance, gives a sig- 
nificant contribution for small radii, but is negligible for 
kp Ro 3> 1. Moreover, the sign of S, which determines 
stability, is essentially fixed by a alone. The relative 
unimportance of the second-order correction is reminis- 
cent of the Strutinsky iheorem^^ for finite fermion sys- 
tems, which states that shell effects are dominated by the 
single-particle contribution in the mean-field potential. 



Stability Diagram 

The stability of a cylindrical nanowire of radius Ro at 
bias V and temperature T is determined by the function 
H(g; Ro, V, T): If E(q) > Vq, then the nanowire is stable 
with respect to small perturbations, and is a (meta)stable 
thermodynamic state. If 3(g) < for any q, then the wire 
is unstable. 

The second term on the r.h.s. of Eq. Ij33(l is posi- 
tive semidefinite, and thus cannot lead to an instability. 
The first term a(q) describes instabilities in two different 
regimes, as in the equilibrium caseiSLS^ (z) The electron- 
shell contribution has deep negative peaks at the thresh- 
olds to open new conducting subbands (c.f. Fig. 01 - 



FIG. 5: (Color online) Stability of cylindrical metal nanowires 
versus Sharvin conductance 124H and bias voltage. Shaded 
(red) areas indicate stability with respect to small perturba- 
tions at T = 0.008 Tf- Solid lines indicate subband thresholds 
for right- and left-moving electrons. 



The surface contribution to a(q) becomes negative for 
qRo < 1, the classical Rayleigh instability. From Eqs. 
(|33|l and (|34p. it is apparent that the most unstable mode 
(if any) within the semiclassical approximation is q = 0, 
except for unphysical radii Ro < 7s/c s (less than one 
atom thick). To illustrate this point, the second deriva- 
tive D(R , V, T) = d 2 E/dq 2 \ =Q is shown in Fig. H Note 
that D > for k F R > 3. 

The stability properties of the system are thus com- 
pletely determined^ by the sign of the stability function 
A(Rq,V,T) ee E(q = 0;R o ,V,T). Figure shows a sta- 
bility diagram in the voltage V and radius Ro plane. The 
x-axis is given in terms of the Sharvin conductance 12411 . 



_ 0.6 



> 
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G s [2eVh] 

FIG. 6: (Color online) Stability of cylindrical metal nanowires 
versus Sharvin conductance Gs between 10 and 50Go and bias 
voltage. Shaded (red) areas indicate stability with respect to 
small perturbations at T = 0.008 Tf- 
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to facilitate the identification of the quantized (linear- 
response) conductance values of the stable nanowires. 
The shaded regions show nanowires that are stable with 
respect to small perturbations, with darker regions rep- 
resenting larger values of A(Ro, V, T). In the figure, the 
solid lines show the subband thresholds for right- and 
left-moving electrons, which are determined by Eq. (1221 . 
At the temperature shown T = 0.0087f, which corre- 
sponds roughly to room temperature, the electron-shell 
effect dominates, leading to instabilities at the subband 
thresholds, and stabilizing the wire in some of the inter- 
vening fingers. 

A stability diagram up to Gs — 50Go is shown in 
Fig. |H1 where the subband thresholds have been omitted 
to avoid clutter. Figures [5] and |S] show that cylindrical 
metal nanowires with certain magic conductance values 
G/G = 1,3,6, 12, 17, 23, 34, 42, . . . remain linearly stable 
at room temperature up to bias voltages eV ~ O.Ief or 
higher. These ballistic conductors can therefore support 
extremely high current densities, of order 10 11 A/cm 2 
by Eq. Q. These are precisely the same magic cylin- 
ders which were previously found to be linearly stable at 
zero bias up to very high temperatures^ Cylinders with 
G/Gq = 8 and 10 are also predicted to be stable at finite 
bias, but not as robust as the neighboring configurations 
with G/Gq = 6 and 12. 

It should be mentioned that in addition to the sta- 
ble cylindrical configurations shown in Figs. |S] and [5] 
nanowires with elliptical cross sections and conductance 
G/Gq = 2,5,9,29,... were also found to be stable at 
zero biasjSi although their finite-bias stability has not 
yet been investigated. 

Perhaps the most startling prediction of Figs.EJandElis 
that there are a number of cylindrical nanowire structures 
which are stable with respect to small perturbations at fi- 
nite bias, but unstable in equilibrium! These metastable 
structures could lead to additional peaks in conductance 
histograms at finite bias, which are not present at low 
bias. It may also be possible to observe switching behav- 
ior between the various stable structures as the voltage 
is varied. 

The results of the above stability analysis should be 
directly relevant for nanowires made of simple monova- 
lent metals, such as alkali metals and, to some extent, 
noble metals. Indeed, the calculated bias dependence 
of the stability of metal nanocylinders with conductance 
G/Gq = 1 and 3, shown in Fig. EI is consistent with ex- 
perimental histograms for gold nanocontacts^ where a 
peak at G « Go was found up to 1.9V at room tempera- 
ture, and a peak at G « 3Go was found up to about 1.5V. 
Similar experimental results have been obtained by sev- 
eral groupsi 3 ' 4 ' 5 ' 6 ' 7 ' 8 Our analysis strongly suggests that 
the remarkable stability properties of gold nanowires at 
finite bias are not a special property of gold, but rather 
a generic feature of metal nanoconductors. 



Nature of the instability 

The nature of the predicted instability of metal 
nanocylinders at finite bias may be illuminated by means 
of a nontrivial identity 22 linking Eqs. I|23[) and (|33|l : 



limH(g;i? 0) V,T) 

g->0 



/ d 2 



1 d \ Sl(Ro,V,T) 



\dR 2 a 



i?o dRo J L 

This implies that the instability corresponds to a 
homogeneous-inhomogeneous transition^ since the r.h.s. 
of Eq. H36|) is proportional to the energetic cost of a 
volume-conserving phase separation into thick and thin 
segments. In the inhomogeneous phase at finite bias, 
the surface corrugation will not be static, but will dif- 
fuse like a defect undergoing electromigration^^^ The 
stable nanocylinders are immune to electromigration, be- 
cause they are translationally invariant and they are so 
thin that they are defect free. Electromigration is possi- 
ble only if a surface-defect is nucleated^ 3 - which becomes 
energetically favorable on the stability boundary. The 
predicted surface instability may thus represent the ulti- 
mate nanoscale limit of electromigration. 



V. CONCLUSIONS 

In this paper, we develop a self-consistent scattering 
approach to the nonequilibrium thermodynamics of open 
mesoscopic systems, and use it to study the cohesion and 
stability of metal nanocylinders under finite bias. In our 
approach, the positive ions are modeled as an incom- 
pressible fluid, and interactions are treated in the Hartree 
approximation, using a quasi-one-dimensional form of 
the Coulomb interaction. This single-band model is ap- 
propriate for simple monovalent metals. It is especially 
suited to alkali metals, but is also appropriate to describe 
quantum shell effects due to the conduction-band s elec- 
trons in noble metals. 

We have utilized a semiclassical treatment of the 
electron-shell structure that plays a crucial role in 
stabilizing metal nanowires at finite bias. Previous 
studies 24,4 — have shown that this semiclassical approach 
accurately describes the energetic cost of long- wavelength 
surface perturbations, which are the limiting factor— in 
the structural stability of long nanowires. Furthermore, 
we have assumed a ballistic nonequilibrium electron dis- 
tribution in the nanowire at finite bias, neglecting in- 
elastic electron-phonon and electron-electron scattering. 
This approximation is valid for wires shorter than the 
inelastic mean-free path. 

We find that the tensile force in a nanowire can be 
modulated by several nano-Newtons when biased by 
a few volts. Such a large effect should be observ- 
able experimentally^ 4 ^*^ although the intrinsic behav- 
ior might be masked by electrostatic forces in the external 
circuit. 

The principal result of this paper is a linear stability 
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analysis of metal nanowires at finite bias, which reveals 
that cylindrical wires with certain magic conductance 
values G/Gq = 1, 3, 6, 12, 17, 23, 34, 42, . . . remain stable 
up to bias voltages eV ~ 0.1 Ep or higher, with the max- 
imum sustainable bias decreasing with increasing radius. 
In particular, wires with G/Go = 1 and 3 are predicted to 
be stable up to eV ~ 0.5 ep- This maximum voltage is 
slightly larger than what is observed experimentally 2,8 
It should, however, be pointed out that stability with 
respect to small perturbations is not a sufficient condi- 
tion for a nanowire to be observed: Metal nanowires are 
metastable structures, and can be observed only if their 
lifetime is sufficiently long on the experimental timescale. 
As a result, the observed maximum sustainable bias is 
likely to be somewhat smaller than that predicted by a 
linear stability analysis. 

A striking prediction of our stability analysis is the 
existence of nanowire structures (e.g. cylinders with con- 
ductance G/Go = 2, 5, 7, 9, 14, 20, . . .) that are only sta- 
ble under an applied bias. This suggests that conduc- 
tance histograms taken at finite voltage might have ad- 
ditional peaks, or even a completely different set of peaks, 
compared to zero- voltage histograms. It may also be pos- 
sible to observe switching between different stable struc- 
tures as a function of voltage. 

Metal nanowires with elliptical cross sections and con- 
ductance G/Go = 2, 5, 9, 29, . . . are also predicted to be 
stable at zero biasi^i Although some of the conductance 
values of the elliptical wires coincide with those of cylin- 
drical wires predicted to be stable only at finite bias, it 
should be possible to distinguish these geometries experi- 
mentally due to the different kinetic pathways involved in 
their formation, and the very different bias dependence 
of their stability. 

Finally, we point out that the predicted instability of 
metal nanowires at finite bias may represent the ultimate 
nanoscale limit of electromigration, due to the current- 
induced nucleation of surface modulation in an otherwise 
perfect, translationally invariant nanowire. 
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Local Density of States gr{z,E) 

In order to include the temperature in the semiclas- 
sical formalism, we use a convoluted density of states 
g T (z,E) = J(-df (E-E')/dE)g{z,E')dE', where 
fo(E) = [l + cxp(f3E)]^ 1 . Thermodynamic quantities are 
then obtained through their zero-temperature expression 
with the density of states g(E) replaced by gr(E). The 
temperature dependence of the average part g~T(z,E), 
proportional to (1 + 0{T/T F ) 2 ) « 1, (Tp is the Fermi 
temperature), is negligible, while the temperature depen- 
dence of the fluctuating part 5gr(z, E) is included in the 
damping factor a vw (T) [see Eq. (TH?1) ]. In the following, 
we set the factor a vw (T) equal to its unperturbed value 
at the Fermi energy sp, since the variation of this fac- 
tor with the perturbation or with energy does not give an 
important contribution. We also drop the subscript T for 
the convoluted density of states to simplify the notation. 

For a perturbed cylinder, the average part of the den- 
sity of states, Eq. iJTSJ, can be expanded to second order 
in the small parameter A as 

,(0)^ A ^(D +A 2^(2) j 



g(z,E)~gl 



(Al) 



where 



E 

E 

Ep 



kpRp 

hp 
4e F 



i(2) _ 



kp Rq 

2-kef 
( kpRp 

kpRo 
6itef 

2tt£ F \ Ep y ' 



6lT£l 



5R{z) 



E ' 



£ -§SR"(z), 



kpRp 
8e f 



SR' 2 (z) 



6ttef 



] j^5R(z)5R"(z) 



and the prime denotes differentiation with respect to z. 

Similarly, the fluctuating part of the density of states 
for a small deformation of a cylinder, Eq. i|19|) , can be cal- 
culated using semiclassical perturbation theory^SLSi and 
is found to be 

6g(z,E)~6gW(E)+\5gW(z,E) + \ 2 6gW(z,E), (A2) 

with 



Sg 



Sg 



(i) 



(o) _ 



hp 



2-kef 



■ cos 9 VW (E), 



fop QvwiT^ f vw L vw c-jj/ \ 



2tte 



APPENDIX A: EXPANSION OF THE 
NON-EQUILIBRIUM FREE ENERGY 

In this appendix, we present more details of the deriva- 
tion of Eqs. ptj )l and 



Sg 



(2) 



x [cos 9 VW (E) - k E L vw sin9 vw (E)] 



kp 

2-KEf, 



^ ^ a>vw 



(T).fvwL 2 w 
2v 2 R 2 



k E SR 2 (z) 



x [2 sin 9 VW (E) + k E L vw cos 9 VW (E)] 
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where, once more, the sum is over all classical peri- 
odic orbits (v,w) in a disk billiard of radius Rq, the 
factor f vw = 1 for v = 2w and 2 otherwise accounts 
for the invariance under time-reversal symmetry of some 
orbits, L vw — 2vRq sin(7rw/i>) is the length of peri- 
odic orbit (v,w), 9 VW (E) = kpL vw ^J E / ep — 3i>7r/2, and 
a vw (T) = t vw / smh.T vw , with t vw = irk F L vw T /2T F , is a 
temperature dependent damping factor. 

To shorten subsequent equations, we define the func- 
tions of energy (E) by writing the first- and second- 
order contributions to the local density of states as 

9 (1 \z,E)=g[ 1 \E)6R(z)+g { 2 1 \E)SR"(z), (A3) 
g^ ( Z ,E) = gW (E)SR 2 (z) + (E)SR' 2 (z) 

(A4) 



+ g {2) {E)5R{z)5R"{z). 



The total density of states per unit length of a cylindrical 
wire is gW(E) = g^(E) + Sg^(E). 



Bare charge imbalances Spo(z) and Spo(z) 

Now substituting the expansion of the local density of 
states into Eq. I|26|l . one gets the expansion of Spo(z) as 



with 



5Pq\z) = Q^8R(z) + q^ ] 8R"{z), 



where the coefficients g^' and g^ are defined as 



.3/2 3/2 



g 0) = k 3 F Ro I £+ +e- 



3tt 



-3/2 



£ F I 2 



(i) k F Ro ( e 1 / 2 + s 1 / 2 \ 
A = — I — 2 | , 



6tt 



-1/2 



(A7) 



with e± = p,± — Uq, and 

(1) kp£± fvwdvw(T) 



5g\ 



ttef * — ' V 2 Rq 

vw 



X i^Lyw cos VW (e±) + k± sin 9 



where k± = k F y/s±/sp, while 

Sp {2) (z) = g (2) SR 2 (z) + g (2) 8R' 2 {z) + g {2) SR(z)SR" \z), 



where coefficients g^ , gp and g^ are 



9? 



(2) 

9 2 



, 3 / 3/2 3/2 
tip I £ i — i £■ 



67T 



-3/2 
~-F 



kpRo f e + + e_ 



16 



^ 2) + ^L 2) 



(i) 



and 



r (2) _ k F e± fvw a vw(T) T 

°9± — ~ „ / , T7v? ^vw 

ZTV £p ^— ' 



x \L VW s\nQ vw {e±) - 2k± cos0 VU) (£±) 
Similarly, Eq. (|29|l for <5po(z) is expanded as 
6p ( Z ,E)~5pf ) +XSp ( 1] +\ 2 5p (2 \ 



(A8) 



with 



(o) _ U 



>)( e+ )+/)( e _) 



(!) - 



JpW = ^ 2) (5i? 2 (z) + 4 2) ^ ,2 (^) + ^3 2) ^) l 5i?"(z), 



(A5) where gf = gf + \U Q \gf (e + ) + gf> 

Effective Coulomb potential V(z, z') 

The expansion of the Coulomb potential 10 as a series 
in A gives 

V{z,z')~V^{z-z') 

+ \V {1) {z,z')+\ 2 V {2) {z 1 z') 1 (A9) 



Mi 



where 



V z2 + vR 

(1) = 6R(z) + SR(z') dVW(z-z') 

[ ' j 2 dHo 

y( 2 )(z,z') = i([^) + ^')] 2 ^ 

[<Si2(«) -<5i?(^')] 2 \ dV< > («-«') 



For future use, let us define the Fourier transform of 
V (0 \z) as 

V(q) = f dze- iqz V {0) {z). (A10) 
Jo 

Note that ReF(g) > 0. 
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Inverse dielectric function e 1 (z,z') 

We first expand the dielectric function e(z,z'), Eq. 
(j?7j) . as 

e(z,z') = e (0) (z-z') 

+ Ae«(z,z') + A 2 e( 2 >(z,z'), (All) 
where the zeroth-order term is 

e(°)(.) = 5{z) + g(0)(£+) + g( ° )(£ -V (°)(z), (A12) 
the first-order term is 

£ W( V ') = ^( S W( Z ,a(»)(,-z') 

+ gW(z,e a )vW(z,z')), 



a=± 



and the second-order term is 

Ml, v ' 



eV{z,z') = \Y.{v {2) ^ a )V^{z-z') 



a=± 



+ (z, 6a )vW (z, z') + gW (z, e a )V^ (z, z')) • 
Let us define the Fourier transform of (z) as 



(g) = / dze-^e (0) (^) 



(A13) 



Note that Re (i(q)/V(q)j > 0, since both .g (0) (e) > 
and V(q) > 0. Substituting Eq. (|AL1|I into the identity 

dz"e- 1 (z,z")e(z",z')=5{z-z'), 

one can solve order by order for the inverse dielectric 
function e^ 1 (z, z') 

e- 1 (z,z') = e- 1 'W(z-z') 

+ \e- 1 'W(z, z') + A 2 e-^ 2 \z, z'), (AH) 



where the zeroth-order term is 



■ 1>(0) (*) = tE 



1 x-^ e I9Z 



the first-order term is found to be 

e- 1 ^ 1 ) = - J "dz x J 'dzze-^iz - z x ) 

xe^(z 1 ,z 2 )e- 1 ^(z 2 -z'), 



and the second-order term is 

e-^ 2 )(z, z') = -J d Zl J dz 2 e^ (0) (z 2 - z') 

- 1 >(°\z-z 1 )eW(z i ,z 2 ) 

+ e- 1 -^(z-z 1 )e^(z ll z 2 ) 

Screened Potential V(z,z') 



Substituting the expansions (| A9|) and (|A14|) of V and 
er 1 into Eq. ijHH , one gets an expansion of the screened 
potential as 

V"(js,jO = V i0 \z-z') 

+ XV il) (z 1 z') + \ 2 V (2 \z,z') (A15) 
where the zeroth-order term is 

V {Q \z) = f dz! V^iz-z^e-^Hz!) = yJ2 tQe** 



L ^ e(q) 



the first-order term is 
V^\z,z') = \dz x 



V^{z-z 1 )e- 1 ^ 1 \z 1 ,z') 

+ V^(z-z 1 )e- 1 ' ( -°\z 1 -z') 



while the second-order term is 



V (2) (z,z') = dzi 



V^iz-z^e-^izuz') 
+ vV>(z-zi)e- 1 'W(z 1 ,z') 
+ V^(z~z 1 )e- 1 ^(z 1 ~z') 



Grand canonical potential Q[V, R(z), U(z)] 

Using Eqs. fKJ). l|A2j). (|X5|) . (|X8l) . and PUHl for 

g(z,E), Spo, 5po, and V, we are now ready to expand 
f2, starting by rewriting Eq. (|30[1 as 

n~n [{v±},R(z),u }-N + u 

- nxK/ii}, c/ ] + « 2 [{m±}, 0b] 

+ fi 3 [{^ ± },i?(z),C/ ], (A16) 

where 



Qi = Uq I dz Spo(z) 
~U 



U L 



^2g^(e a )+Xg^b(0) 



a=± 



a'E^ + c^-sSV] i%)i 
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O, = 



(Uo) 



(U f 



^2 J dzg(z,e a ) 



a=± 



a=± 



9^(ea)+XgP(e a )b(0) 



^[^ 2) ( £Q ) + [ 5 ( 2) ( £Q )- 5 f( £Q )] g 2 ] |6(g)| 



and ^3 = i J dzdz'5po(z)V(z, z')Spo(z) can be written 



+ A [{/x±}, iZo, C/ ] + A 2 ^ J [{ M± }, i? , C/o] 
The zeroth order term in the expansion of ^3 is 



2o(2)f 



(0) 



- / dzdz' Sp ( ° ] [z)V<® (z - z')8p { ^ {z') 



a=± 

the first-order term is 



Oi lj = x / dzdz' 



1 



4 O) (z)T>(%,z0^1V) 

+ 2^ 1 t^ (0 ^-^)4V) 

= LU Q i 1 1) b(0) + O(L ), 



and the second-order contribution is 



(2) 



dzdz' ^ 0) {z)V (2) (z, z')Sp ( ° } (z') 

+ 26ti\z)V^(z-z>)6p ( °\z>) 



L v V(q) 
2 ^ 



(1) (1) 2 
01-02 9 



>( g )| 2 +O(L ). 



Adding up all the contributions in Eq. (|A16|) . and drop- 
ping contributions of order L°, one gets Eqs. (|S2l and 

The above calculations also show that Eq. I|30(l can be 
rewritten as 



n = n a [{ f i ± },R(z),u a ]-N + u 

+ I [ dzdz' Sp ( ( p (z)VW (z - z')Sp^ (z') 



+ O{L ). (A17) 
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